Stable trimorphic coexistence in a lattice model of 
spatial competition with two site types 



Ilmari Karonen* 

Department of Mathematics and Statistics, University of Helsinki, Finland 



Abstract 

I examine the effect of exogenous spatial heterogeneity on the coexistence 
of competing species using a simple model of non-hierarchical competition 
for site occupancy on a lattice. The sites on the lattice are divided into two 
types representing two different habitats or spatial resources. The model 
features no temporal variability, hierarchical competition, type-dependent 
interactions or other features traditionally known to support more compet- 
ing species than there are resources. Nonetheless, stable coexistence of two 
habitat specialists and a generalist is observed in this model for a range of 
parameter values. In the spatially implicit mean field approximation of the 
model, such coexistence is shown to be impossible, demonstrating that it 
indeed arises from the explicit spatial structure. 
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1. Introduction 

The coexistence of competing species, and factors promoting and limiting 
it, are of considerable practical and theoretical interest in ecology. A well 
known "rule of thumb" is the principle of competitive exclusion (Gause, 1932; 
Tilman, 1982; Levin, 1970), which states that at most n mutually competing 
species may stably coexist on n available resources. 

Under suitable assumptions, the competitive exclusion principle can be 
proven as a mathematical theorem. However, if these assumptions are vi- 
olated — for example, if the resource abundances may fluctuate over time, 
either due to external resources or simply because the ecological dynamics 
tend to a cyclic or chaotic attractor — it may no longer hold (although a 
related concept, the essential dimensionality of the environment (Dieckmann 
and Metz, 2006; Metz et al., 2008), may still be applied to such systems). 

For systems which do satisfy these assumptions, the validity of the com- 
petitive exclusion principle depends fundamentally on just what we count as 
"a resource" (Abrams, 1988). This is not as simple a matter as it sounds. 
Were one to ask a practical-minded ecologist what constitutes a resource, 
they might name examples such as water, sunlight and nutrients for plants, 
or prey species for animals. But in the mathematically exact form of the 
competitive exclusion principle, almost anything may constitute a distinct 
resource: a single prey individual, a square meter of land, a speciflc com- 
bination of nutrient concentrations, etc. Thus, even in systems which the 
competitive exclusion principle formally holds, the actual number of poten- 
tially coexisting competitors may be greater than one would expect by naively 
under counting the resources. 

The effect of spatial structure on the maximum diversity a system can 
support is, in particular, often neglected. For example, systems consisting 
of several distinct types of habitats are often modelled by assuming that 
each habitat constitutes a homogeneous patch, within which the populations 
are well mixed. If competition between individuals is for suitable living 
space within these habitats, space in each habitat then becomes a single 
resource, and thus one would expect (and will, given these assumptions, 
mathematically discover) that at most n competitors may stably persist in 
n distinct habitats. 

In reality, however, even if habitats are homogeneous, they are certainly 
not usually well mixed. Thus, individuals living near habitat boundaries 
will, over time, experience a different environment than those living in the 
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interior of habitat patches. (Even if the individuals themselves don't move 
or interact with anything outside their local site, their offspring must still 
disperse and may end up in a different habitat.) This can create additional 
niches near habitat boundaries in which additional competing species might 
be able to coexist. 

To demonstrate this, I present in this paper a simple spatially explicit 
toy model of site occupancy competition, which supports stable coexistence 
of three strains - two specialists and one generalist - on two spatially seg- 
regated habitats on a lattice of sites. This model contains no other features 
known to promote coexistence, such as internally or externally generated 
temporal fluctuation (Hsu et al., 1978; Armstrong and McGehee, 1980), hi- 
erarchical site competition (Adler and Mosquera, 2000; Tilman, 1994), direct 
strain-dependent competition terms (Murrell and Law, 2003) or cooperative 
or other nonlinear interactions between individuals. Rather, the coexistence 
arises merely from the presence of habitat boundaries combined with passive 
distance-limited dispersal, which causes specialist strains to be locally mal- 
adapted near these boundaries and thereby allows more generalist strains to 
persist there. 

As far as I know, this particular mechanism of coexistence has not been 
previously studied in individual-based models. A similar coexistence mech- 
anism was very recently described in a reaction-diffusion model and in a 
ID stepping stone model by Debarre and Lenormand (2011), who termed it 
"habitat boundary polymorphism". The results in this paper parallel theirs, 
and confirm that this mechanism is robust with respect to the details of the 
model, provided that the basic features of habitat heterogeneity and passive 
distance-limited dispersal are present. 

A model almost identical to mine was studied by Lanchier and Neuhauser 
(2006), who showed analytically that it could support stable dimorphic co- 
existence, either of a generalist and a specialist strain or of two different 
specialist strains. My model differs from theirs only in that they restrict the 
habitat configuration to the special case of an infinite regular checkerboard 
pattern consisting of n-by-n squares of each habitat type.^ However, while 
Lanchier and Neuhauser did briefly remark that "the generalists persist for 
a very long time along the boundaries . . . where the density of specialists 



^Lanchier and Neuhauser (2006) also consider a somewhat different set of dispersal 
kernels, but both theirs and mine include the basic case of strict nearest-neighbor dispersal. 
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is low" in numerical simulations, they do not seem to have investigated this 
possibility of trimorphic coexistence in their model further. Similarly, Snyder 
and Chesson (2003) define a model quite similar to mine (although in discrete 
time and one spatial dimension), and observe the enhancing effect of stable 
spatial heterogeneity and local dispersal on coexistence of competitors, but 
also restrict their analysis to two competitors. 

2. Model definition 

I model a population of haploid, asexually reproducing sessile organisms 
with distance-limited offspring dispersal. The model I'll define below belongs 
to the class of lattice contact processes (Harris, 1974; Neuhauser, 1992), in 
which the environment is taken to consist of a lattice of discrete sites, and 
interactions are (mainly) between nearest neighbor sites. 

Let L be a regular two-dimensional lattice of sites (e.g. L = Z^, although 
for numerical simulations a finite lattice must obviously be used). To each 
site I assign a random, fixed habitat class ("A" or "B"), such that both classes 
of sites are present in L in equal numbers. (I will describe the way in which 
the habitat classes are assigned in more detail below.) 

Each site in L may, at a given time, be either vacant or occupied by an 
individual belonging to one of three strains: an A-specialist (a), a B-specialist 
(6) or a generalist (g). The two specialist strains can only occupy sites in 
their respective habitat class, while the generalist strain may occupy any site. 

Except for their different habitat adaptation, the strains are completely 
identical: all individuals die with rate fi and produce offspring with rate (/). 
Offspring are randomly dispersed to the eight nearest sites surrounding the 
parent individual's site (or possibly, with probability e, to a randomly chosen 
site in L), and will become new individuals if the site they land in is vacant 
and of a suitable habitat class. However, the generalists pay a cost for their 
ability to live in either habitat: their offspring survive only with probability 
Pg < 1.' 

The time evolution of the entire lattice L can thus be considered as a 
continuous-time Markov process, whose state at time t is a function rjt : L ^ 
{0, a, b, g} mapping sites in L to their occupancy state (with the state 



^ Equivalent ly, I could' ve scaled the offspring production rate of the generalists to 
Pg(j). However, the definition I've chosen allows a straightforward generalization to semi- 
specialist strains with different (non-zero) survival rates in different habitats. 
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denoting a vacant site). The local transition rates of a site x are then 
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r(a ^ 0) = r(b ^ 0) = r(g -> 0) = /i, 

where r{s s^) is the rate at which a site in state s changes to state 5', 
is the habitat class of site x and Ex is the set of sites adjacent to x. 

The behavior of the model is determined by the habitat configuration 
along with two parameters: the scaled baseline fecundity rate \ — (j)/ 11 {or 
equivalently, the scaled mortality rate = stnd the generalist survival 
rate p^. The baseline fecundity affects the equilibrium population density of 
both the specialists and the generalist strain: at low A all strains die out, 
while at very high A almost all sites are occupied at any time. 

The parameter < ]9g < 1 determines the penalty which the generalist 
must pay for its ability to exploit both site types, and is (together with the 
habitat configuration) crucial in determining the outcome of the model. If 
Pg is too low, the generalist strain will not be viable, or, even if viable on its 
own, may lose in competition to the specialists. Conversely, if is close to 
1, the specialist strains gain little or no advantage over the generalist from 
their specialization, while paying a considerable price in being able to live in 
only one habitat, and can thus be expected to lose to the generalist. 



3. Landscape generation 

An issue so far overlooked in the model definition above is the way in 
which the lattice sites are assigned to their habitat classes. The simplest 
way to do so, of course, is to simply assign each site independently to either 
habitat with equal probability. This produces a lattice with no correlations 
between the habitat classes of different sites. 

However, real environmental features are often correlated, and it would 
be desirable to consider the effects of such correlations on the behavior of the 
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model. To first order, such correlations can be characterized by the pairwise 
correlation probability k = Pt[Hx = Hy\y G Ex]^ i.e. the probability that 
two randomly chosen adjacent sites have the same habitat type. With an 
equal number of sites in each habitat, the pairwise correlation probability on 
a completely random, uncorrelated lattice is A; = |, while lattices with | < 
k < 1 are positively correlated and those with < < | are anticorrelated.^ 

To generate random habitat class distributions with a given pairwise cor- 
relation probability for numerical simulations, I start by randomly choosing 
half of the sites and assigning them to habitat A and the rest to habitat B. I 
then apply an iterative annealing process, which reassigns randomly chosen 
sites to new habitat classes with suitable weighted probabilities, until the 
desired value of k is reached. 

There are many possible variations of the general annealing scheme I've 
used. The basic idea in all of them is to change the habitat classes of ran- 
domly chosen sites if this would change k in the desired direction, while 
occasionally also allowing changes in the other direction so that the process 
doesn't get stuck at a local minimum or maximum. 

The particular annealing algorithm I've used to generate habitat config- 
urations for the simulations in this paper^ chooses random adjacent pairs of 
sites, and swaps their habitat classes with probability 

- 

^~ d^ + {l- d)^ 

if ^current < ^target, 01" with probability 1—p otherwise, where d is the fraction 
of sites adjacent to the chosen pair which belong to the opposite habitat than 
their neighbor in the chosen pair. 

The exponent 7 is a free parameter which controls the "temperature" of 
the process. When 7 = 1, the probability of exchanging the habitat classes 
of a chosen site pair is a linear function of d. This tends to produce fairly 
slow convergence and rough, jagged cluster boundaries. At high values of 7, p 
approaches a step function, producing faster initial convergence and smoother 
cluster boundaries, but also increasing the risk of the process getting stuck 
at a local maximum or minimum of k. 



^On a square lattice where each site is adjacent to its 8 nearest neighbors, the smallest 
achievable value of /c is \. 

^The interactive Java applets from which the snapshots in figure 6 are taken use a 
different annealing algorithm. 
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4. Mean field approximation 

Classical ecological theory predicts that the only possible outcome of three 
distinct strains competing for the occupancy of two habitats should be the 
eventual extinction of one or more of the strains, except at specific degenerate 
choices of the model parameters where neutral coexistence may occur. I will 
demonstrate below that this prediction indeed holds if the populations are 
assumed to be well mixed, either globally or within each habitat. However, I 
shall also show that, in the full model with explicit spatial structure, a region 
of stable trimorphic coexistence does exist for intermediate values of Pg. 

Assuming that all offspring disperse uniformly over the entire lattice, i.e. 
that € = 1 in 1, the transition rates of each site are fully described by the 
mean population densities na, nb and of the diflFerent strains, where 



for each strain s. Further letting the lattice size \L\ tend to infinity, one 
obtains a simple system of ordinary differential equations describing the 
time evolution of these mean population densities — the so called mean field 
approximation — which may be solved analytically. Such an approximation 
of an essentially equivalent model was presented by Lanchier and Neuhauser 
(2006), who showed that trimorphic coexistence was only possible for degen- 
erate choices of parameter values. 

However, simply assuming all dispersal to be global completely neglects 
not only the detailed spatial habitat structure, but even the pairwise corre- 
lation parameter k. A better approximation, similar to the "improved mean 
field approximation" of Hiebeler and Morin (2007), is obtained by assuming 
well mixing only within each habitat. The resulting approximation can be 
interpreted as a model of a population inhabiting two well-mixed habitat 
patches, with a fraction k of all offspring remaining within their parent's 
patch and the rest dispersing to the other patch. This two-patch approxima- 
tion takes into account the habitat correlation parameter k but still retains 
the analytical tractability of the mean field approximation. (For = |, the 
two approximations are exactly equivalent.) Below, I will analyze this ap- 
proximation of the model defined above, and show that it also only supports 
non-degenerate coexistence of at most two strains. 

Assume that the occupancy states of the lattice sites are independent, 
and that the probability of a site being occupied by a given strain s is equal 
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Figure 1: The outcome of the model as predicted by the two-patch approx- 
imation for different values of the normalized mortality rate ji/cf) and the 
generalist survival probability Pg. The diagonal line dit fi/cj) = Pg and the 
vertical line at = mark the boundaries at which the generalist and 
specialist strains go extinct even in the absence of competitors. In the re- 
gions marked "generalist wins" and "specialists win", all strains can survive in 
the absence of competitors, but from any initial state with all three strains 
present, the system converges to either a monomorphic generalist-only state 
or a dimorphic specialist-only state. Along the white line at Pg = fc, the three 
strains can coexist neutrally. 
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to the mean population density Us^h of that type in the site's habitat H. 
Then, in the hmit as |L| ^ oo, the time evolution of the population densities 
can be written as 

Ps,A^A(l){kns^A + (1 - k)ns^B) - H^Us^a, 
Ps,BVB(l){kns^B + (1 - k)ns^A) - l^ris^B^ 

for s G {a, b,g}, where va and vb are the vacant site densities in the two 
habitats and Ps^A and Ps^B are the probabilities of an offspring of type s 
surviving in the respective habitats: 

Pa,A =Ph,B = 1, 
Ph,A = Pa,B = 0, 
PgA = Pg,B = Pg. 

Equivalently, this system may be written in matrix form as 

d 

-—ris — Maris 
dt 

for s G {a, b,g}, where fig — [^s,A5^s,b]^ and 

^ ^ r (t)kps^AVA - /i 0(1 - k)ps^AVA 

^ |_ 0(1 - k)ps^BVB (t>kps,BVB - _ * 

If this system has a nontrivial interior equilibrium n^, this necessarily implies 
that Msfis = [0, 0]^ 7^ n^, and therefore that Mg must be singular, and thus 
have a zero determinant, for each strain s present in the population. Writing 
out the determinant as 

\Ms\ = (t)^{2k - 1)Ps,aVaPs,bVb 

- 4^fik {ps,aVa + Ps,bVb) + /i^ = 0, 

yields, for each 5, an equation containing the same two unknown variables: 
Va and vb- As the coefficients Ps^A and Ps^B will, in general, be different for 
each strain 5, one can see that, except for degenerate choices of the parameter 
values, no solution will exist for more than two strains. 

More specifically, we can see that, in the absence of specialists, small 
generalist populations can grow in the two-patch approximation of this spe- 
cific model if and only if pgcj) > /i, and conversely that small populations of 
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either specialist strain can grow in the absence of the generahst if and only 
if kcj) > ji. Outside these regions, shown in figure 1, the respective strains 
are not viable and will always die out (as the per capita growth rates in this 
model are always maximized at vanishing population densities). 

Within the region where all strains are viable, the approximated model 
has (in general) two possibly stable equilibria: one where only the generalist 
is present, and one where the generalist is absent and both specialists present. 
(Any equilibria with only one specialist present are obviously unstable against 
invasion by the other specialist.) The former is stable if and only if > fc, 
while the latter is stable if and only if < k. Only at exactly Pg = shown 
as the white line in figure 1, both of the equilibria are neutrally stable, and 
are connected by a line of neutrally stable equilibria along which neutral 
coexistence can occur. 

5. Simulation results 

Studying the dynamics of the full, unapproximated model requires numer- 
ical simulations. As such simulations tend to be computationally intensive, 
I have carried them out using custom, optimized programs written in the C 
programming language.^ The simulation code used for this paper includes 
two variants of the coupling-based simulation algorithm described in (Karo- 
nen, in prep.), one using an occupancy list for low population densities, and 
another using a vacancy list for high densities, with the outer simulation 
loop periodically checking the population density and switching to the vari- 
ant with the higher mean time step per iteration. I have also ported the basic 
simulation code (without the coupling technique) to Java for demonstration 
purposes using interactive applets. 

All simulation runs for this paper were done on a square 256 x 256 lat- 
tice with 8 neighbors per site and with the edges wrapping around to the 
opposite sides. In all simulation runs, time has been scaled so that the per 
capita mortality rate /x = 1; in effect, I measure time in mean individual 
lifetimes. Habitat configurations were generated using an annealing method 
as described in section 3. The "fiea" pseudorandom number generator (Jenk- 
ins, 2007) was used to produce random numbers, although I also carried out 
tests using other random number generators to check that the results did not 
depend on such details. 



^Source code available from author under an open source license. 
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Figure 2: Plots of equilibrium specialist and generalist densities at = 0.5 
and = 8/x as functions of Pg obtained from 20 numerical simulation runs. 
See text for details. The graph on the right has been plotted on a logarithmic 
scale and zoomed in to better show the coexistence region. 

Figure 2 shows the equilibrium densities of generalists and specialists ob- 
served in repeated individual-based simulations of the model on a 256 x 256 
lattice with reproduction to 8 nearest neighbors and wrapped edges with an 
uncorrelated habitat distribution and a moderate value of (/)///. Each simu- 
lation run was started from a random habitat configuration and a random 
initial state, with half the sites occupied by generalists and half by special- 
ists. Populations were allowed to equilibrate for 50000//X time units, after 
which population densities were averaged over another 50000//X time units. 
The specialist occupancy fractions are summed over both specialist strains. 

Contrary to the predictions from the mean field approximation, a non- 
degenerate region of the parameter space where all three strains stably coexist 
can be seen in figure 2. This region is displayed more extensively in figures 
3 and 4, which plot the observed region of coexistence against /ll/c/) and Pg 
for the various habitat configurations (anticorrelated, uncorrelated and two 
positively correlated patterns) shown in figure 5. Figure 3 shows results for 
pure nearest-neighbor dispersal (6 = 0), while in figure 4, 1% of all offspring 
were permitted to disperse uniformly over the whole lattice (e = 0.01). 

For figures 3 and 4, each simulation run was started from a random initial 
population on the fixed pregenerated habitat landscapes shown in figure 5. 
Populations were allowed to equilibrate for 20000//X time units, after which 
population densities were averaged over 2000//X time units. The red areas 
show where only the specialist strains survived, while in the blue regions only 
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Figure 3: Results of numerical simulations on a 256 x 256 site lattice as 
functions oi ji/cj) and Pg on the habitat patterns from figure 5 with no global 
dispersal (6 = 0). See text for details. Compare with figure 4 and with the 
mean field predictions from figure 1. 
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Figure 4: Results of numerical simulations on a 256 x 256 site lattice as 
functions oi ji/cj) and on the habitat patterns from figure 5 with occasional 
global dispersal (e = 0.01). See text for details. The marks on figures 4b and 
4d show the parameter values used for the invasion simulations in figure 7. 
Compare with figure 3 and with the mean field predictions from figure 1. 
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(a) k = 0.3 



(b) k = 0.5 




(c) k = 0.75, 7 = 10 (d) k = 0.75, 7 = 3 



Figure 5: The habitat configurations used for the simulations shown in figures 
3 and 4. White squares correspond to habitat A, black squares to habitat B. 
The insets show a 32 x 32 region magnified. All lattices were generated from 
the same random initial state (which is nearly identical to lattice 5b; only a 
very small amount of annealing was needed to make k exactly 0.5) using the 
annealing method described in section 3. Lattices 5c and 5d have the same 
pairwise correlation number k = 0.75, but the different annealing parameters 
used to generate them lead to visibly different higher order correlations and 
to corresponding differences in population dynamics. 
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the generalist strain remained. The hghter shaded area between them shows 
the part of the parameter space where both speciahsts and generahsts sur- 
vived, with the color gradient shown above the figures indicating the relative 
average population densities. 

In both figures 3 and 4, the region of stable coexistence can be seen as a 
more or less wedge-shaped area starting from the point where the viability 
boundaries of the strains intersect, which is where the two-patch approxi- 
mation would predict a line of neutral coexistence (see figure 1). The main 
diflFerence between the figures can be seen in the lower left side of the coex- 
istence region: with global dispersal, the lower boundary of the coexistence 
region is quite sharp, whereas with no global dispersal and high baseline fe- 
cundity the generalist can often survive in small numbers (shown as a 
light orange hue in the plots) even where the specialist dominates. 

This happens simply because the high fecundity allows even small isolated 
population clusters to survive for a long time, and because the strictly local 
dispersal prevents the specialists from recolonizing isolated habitat patches. 
If such a habitat patch happens to end up with no specialist individuals (ei- 
ther because all happen to die out, or because none were present initially), 
the patch can be colonized by generahsts, which are then safe from com- 
petition there. Allowing a fraction of offspring to disperse globally lets the 
specialist strains recolonize such patches, eliminating this effect. 

It can also be seen that the addition of global dispersal generally reduces 
the width of the coexistence region somewhat, although (except for the iso- 
lated patch effect noted above) the reduction is not yet very large for e = 0.01. 
This is to be expected, given that at e = 1 the coexistence region reduces to 
a line, as shown by the mean field (and two-patch) approximation above. 

The results shown in figures 2, 3 and 4 were calculated using a simulation 
technique based on monotone coupling (Karonen, in prep.), which allows the 
system to be efficiently simulated for all values of the parameter in parallel. 
Each line in figure 2 and each vertical stripe (out of 1024 per plot) in figures 
3 and 4 corresponds to one simulation run. Because the simulation technique 
used causes the effects of random demographic fiuctuations on populations 
with different values of Pg within the same run to be correlated, the results 
show stronger correlations within each run than between runs, which should 
be kept in mind when interpreting the figures. 

Figure 6 contain snapshots of simulations run on lattices with different 
site type patterns. It can be seen that, when sites of the same type are 
strongly clustered, large contiguous clusters are dominated by the respective 
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(a) k = 0.5, j9g = 0.43 (b) k = 0.75, = 0.7 (c) k = 0.95, = 0.93 



Figure 6: Snapshots of the population equihbrium after a few hundred mean 
hfetimes for various values of fc, with near the middle of the coexistence 
region. Snapshots were taken from simulations run on a 128 x 128 lattice 
with = 4/x and e = 0.001. The red and green sites are occupied by a 
and b specialists respectively, the blue sites are occupied by generalists and 
the gray sites are vacant. For the blue and gray colors, darker shades are 
used for habitat A and lighter shades for habitat B. These snapshots have 
been taken from interactive Java applets available at http://vyznev.net/ 
ca/coex2env/. 
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specialist strain, while the generalist strain is able to persist in areas near 
cluster boundaries and in isolated minor clusters too small to support a stable 
specialist population. 

On the other hand, when site types are uncorrelated, a different pattern 
is observed. Such lattices contain no large contiguous clusters that could be 
dominated by one specialist strain; instead, the two specialist strains tend to 
occur together in regions where the random distribution of site types happens 
to favor one or both of them. Through competition with the generalist strain, 
the two specialist strains indirectly support one another, even though there 
is no direct interaction between them. 

6. Mutual invasibility 

To a skeptical mind, the results presented above may yet leave some 
doubt about whether the apparent coexistence observed in these simulations 
is indeed genuinely stable, or merely an artifact of slow convergence and 
insufficient simulation time. After all, if the simulation was run for long 
enough on a finite lattice, eventually one of the strains (and eventually all 
of them) would almost surely go extinct simply due to demographic stochas- 
ticity. Thus, it may not even be entirely clear what "stable coexistence on a 
finite lattice" should actually mean. 

On an infinite lattice, as in Lanchier and Neuhauser (2006), a set of 
strains may be said to coexist stably if they can all persist indefinitely long 
with non-zero probability. By this definition, no stable coexistence (or even 
just existence) is possible in a finite system. However, since there do exist 
known results that relate the scaling of the expected time to extinction on 
a finite lattice, as a function of lattice size, to the limiting behavior of the 
model on an infinite lattice, one might be inclined to try and use such scaling 
laws to extrapolate stability from small lattices to the infinite limit. 

This is not the approach I have taken. After all, real habitats and popula- 
tions, like the simulations employed in this paper, are finite — in appealing to 
a definition of coexistence that only works for infinite populations, one ends 
up obscuring the fact that, in reality, if a population of tens of thousands of 
individuals persisting over equally many generations is not considered stable, 
it's hard to say what should be. 

Rather, I wish to demonstrate the stability of the trimorphic coexistence 
in my model in a more direct manner, by showing that it exhibits mutual in- 
vasibility. That is to say, if a small number of individuals of any of the three 
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strains are introduced into a stable population consisting solely of the other 
two strains, the introduced strain will, with positive probability, survive and 
grow in number up to its equilibrium density in the trimorphic equilibrium, 
with the initial part of the growth curve appearing approximately exponen- 
tial. 

Figure 7 shows some simulations demonstrating mutual invasibility at 
six points within the coexistence region of the parameter space. With 100 
initial invader individuals, the invader strain survived and established itself 
in all runs carried out — evidently the invasion probability of a single invader 
exceeds 1/100 at all the sampled parameter values. The population density 
of the invading strain over time shows a distinctive sigmoid shape, with the 
initial growth being approximately exponential. A notable feature visible in 
the plots is that invasion by a specialist strain also increases the population 
of the other specialist strain already present; this happens because both 
specialists are in competition with the generalist strain. The relatively high 
variance seen in some of the plots during the growth phase is due to initial 
demographic stochasticity affecting the time until exponential growth sets 
in; once properly started, the shape of the growth curve is very similar in all 
runs. 

Were the coexistence observed in this model merely neutral, the popu- 
lation density of a newly introduced strain would be as likely to decrease 
as to increase as the result of stochastic fluctuations. The fact that, at the 
sampled parameter values, small populations of each strain instead show a 
clear increasing trend conflrms that this model exhibits true, non-neutral 
coexistence. 

7. Discussion 

In this paper I've demonstrated, using a simple toy model of competitive 
population dynamics on a lattice, that spatial heterogeneity is one of the 
mechanisms by which the competitive exclusion principle can be violated. 
The fact that this cannot occur in well-mixed populations shows that popu- 
lation viscosity and explicit spatial structure are essential to this mechanism. 

Had the model included more than two habitat types, temporal variation, 
hierarchical competition or nonlinear interactions between individuals, the 
coexistence of multiple strains would not have been at all surprising. Yet 
it has none of these, and can still support more than two strains in stable 
coexistence. All that allows such coexistence to persist in this model is 
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Figure 7: Simulations showing mutual invasibility at the parameter values 
marked in figures 4b and 4d; = 8/x for all simulations. A dimorphic popula- 
tion is allowed to equilibrate for 400//X time units, after which 100 individuals 
(i.e. initial density ^ 0.0015) of the third strain are placed randomly on the 
lattice. In simulations of specialist invasion, the invading strain is without 
loss of generality taken to be b. The population densities shown in the graphs 
are smoothed over 5//x time units and averaged over 10 independent simu- 
lation runs; invasion does succeed in all runs. The red, green and blue lines 
show the average a, b and g strain densities respectively, while the solid, 
dashed and dotted lines correspond to different values of Pg within the coex- 
istence region. (In populations with both specialist strains present, the red 
and green lines tend to overlap.) The thin dotted lines are drawn one sample 
standard deviation above and below the corresponding average line. 
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the combination of environmental variation, persistent spatial structure and 
distance-limited dispersal; eliminating any of these reduces the model to one 
capable of supporting no more strains than would be predicted by a naive 
application of the competitive exclusion principle. 

Real organisms do not usually live in a completely homogeneous envi- 
ronment, nor do most of them disperse uniformly over their entire habitat. 
It is obvious and commonly acknowledged that environmental variation can 
increase diversity, yet the fact that, when combined with distance-limited 
dispersal, this increase can be more than linear seems to have attracted little 
attention. Yet the ubiquity of habitat edges and fragmented landscapes in 
nature suggests that it should be possible to find examples of this type of 
coexistence in nature, and indeed that such "edge effects" may contribute to 
the generation and maintenance of ecological diversity in many, if not most, 
ecosystems. 

I find, however, that in many ways this work has raised more questions 
than it has answered. For example, an obvious question would be whether 
the model allows the stable coexistence of more than three strains. An- 
other natural question is whether the coexistence of three or more strains in 
this type of model can also be evolutionarily stable, and further, whether it 
might arise from a mono- or dimorphic state through evolutionary branching 
(Geritz et al., 1998; Magori et al., 2005). Based on limited simulation exper- 
iments, the answer to all of these questions appears to be "yes", although the 
conditions still need to be explored more thoroughly. 
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